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Accomplishments 

1.  A  3D  code  is  developed  for  propagating  the  small  scale  return  signals  to  the  far  field. 

2.  Validated  using  analytical  signals. 

3.  Inner  field  scattering  from  a  wind  turbine  and  aircraft  is  computed 

4.  An  interface  is  developed  to  couple  the  above  far  field  code  with  inner  field  code. 

5.  Far  field  scattering  from  a  single  wind  turbine  is  done  by  coupling  the  far  field  code  to  inner 
field  data. 

6.  Far  field  scattering  from  multiple  wind  turbines  with  and  without  aircraft  is  computed 

7.  Multiple  return  signals  due  to  ground  reflection  are  also  computed. 


Introduction 

The  current  STTR  project  involves  developing  a  method,  based  on  using  modern  radar,  enhanced 
with  the  new  technique,  Wave  Confinement  [WC],  to  accurately  locate  low  flying  aircraft 
approaching  areas  such  as  Air  Force  bases.  Conventional  schemes  involve  estimating  the  target 
location  based  on  assumptions  of  the  radar  propagation  path  that  often  neglect  important  effects, 
when  integral  equation  methods  are  used.  These  effects  include  index  of  refraction  variations, 
terrain  reflections  and  "clutter"  from  nearby  objects.  Currently,  the  only  alternative  method — "Ray 
Tracing” — can  take  some  of  these  effects  into  account,  but  it  involves  an  incoherent  collection  of 
"rays",  from  which  it  is  difficult  to  extract  information.  Ray  tracing  computations  can  even  become 
chaotic  near  caustics,  or  neglect  caustics  altogether.  Details  of  the  basic  method  are  included  in  Ref1 
and  some  new  developments  are  briefly  discussed  in  this  report. 

The  consequences  of  windmill  farms  to  radar  systems  presents  a  difficult  challenge,  especially 
when  windmill  farms  are  located  near  Air  Force  bases.  Unfortunately,  the  requirements  for  locating 
the  bases  and  farms  results  in  conflicts,  when  the  bases  in  flat  areas  away  from  populated  areas  and 
wind  farms  on  higher  ground  surrounding  them  (as  is  often  the  case].  The  problem  that  we  are 
addressing  concerns  any  low  flying  aircraft  approaching  a  base  over  such  a  wind  farm. 

Characterization  of  radar  returns  from  a  windmill  farm  would  be  feasible  if  the  individual  windmills 
were  placed  far  apart,  or  were  all  in  air  flowing  at  the  same  speed  and  direction.  However,  the 
farms  often  have  several  hundred  turbines  spaced  only  a  few  diameters  apart,  each  of  which  is 
exposed  to  varying  winds  because  they  are  at  different  locations  in  a  complex  topography.  The 
scattered  signals  will  undergo  refraction,  multiple  reflections  and  diffraction.  In  the  present 
project,  diffraction  is  not  addressed  although  some  work  has  been  done2. 

The  final  outcome  achieved  at  the  end  of  the  phase  I  contract  is  the  demonstration  of  the  capability 
to  compute  a  return  signal  from  a  multiple  wind  turbines,  when  an  aircraft  is  flying  over  as  shown 
in  Figure  1.  This  capability  can  be  used  in  future  to  a]  determine  the  location  of  any  aircraft  present 
and  b]  develop  a  tool  to  model  the  atmosphere,  which  plays  a  major  role  in  accurate  detection. 


Receiver 


Figure  1:  Problem  to  be  solved  during  this  project 


iSteinhoff,  J.  and  Chitta,  S.,  “Solution  of  the  scalar  wave  equation  over  very  long  distances 
using  nonlinear  solitary  waves:  Relation  to  finite  difference  methods,”  Journal  of  Computational 
Physics,  Vol.  231 ,  No.  19,  2012,  pp.  6306-6322. 

2  Steinhoff,  J.  and  Chitta,  S.,  “Long  distance  wave  computation  using  nonlinear  solitary  waves,” 
Journal  of  Computational  and  Applied  Mathematics,  Vol.  234,  No.  6,  2010,  pp.  1826-1833. 


There  are  two  spatial  scales  that  must  be  resolved  for  the  above  wave  propagation  problems: 


a  Large  Scale  -  On  the  scale  where  variations  in  the  medium  and  scattering  surfaces  are  0(L], 
the  wave  packets  appear  as  propagating  and  reflecting  codimension  one  surfaces.  (2-D  for 
3D  solutions  for  thickness,  d<<L.]  These  surfaces  propagate  along  the  local  normals  (for 
isotropic  media],  refract  and  undergo  reflection.  (Only  specular  will  be  considered  here.] 
Thus,  at  this  scale,  the  problem  involves  only  surfaces  either  stationary  (reflecting]  or 
propagating  (waves]  along  local  normals. 

b  Small  Scale  -  The  small  scale  wave  packets,  just  like  the  above  surfaces,  propagate  along 
"rays"  as  in  Hamilton  Jacobi  methods.  The  "waveform"  that  they  carry  is  invariant  in  time3 
except  for  dissipation  and  geometric  changes  in  intensity.  The  latter  can  be  included  as 
multipliers  of  the  wave  forms,  so  that  only  initial  wave  forms  need  to  be  specified  for  each 
ray. 


Methodology 

In  this  report,  we  describe  a  new,  simple,  computationally  efficient  wave  simulation  method  that 
overcomes  most  of  the  problems  of  conventional  methods  for  long  distance  propagation  in 
inhomogeneous  media,  including  multiple  reflections.  In  its  initial  form,  it  involves  the  scalar  wave 
equation  with  non-dispersive  and  non-diffusive  media,  but  these  limitations  can  be  removed  in 
future  versions  with  perturbation  terms.  This  new  method  is  termed  the  "Generalized  Eikonal 
Method"  (GEM], 

The  solution  at  any  distant  receiver  ("distant"  meaning  "many  wavelengths  away"]  is  assumed  to  be 
represented  by  smooth  variables  (away  from  caustics]:  attenuation  factor  ( Ak  ],  propagation  vector 

[Sk  ],  arrival  time  (t].  It  is  assumed  that  the  wave  equation  is  accurately  computed  (with 
conventional  methods]  in  a  small  region  surrounding  the  source/target,  with  dimensions 
comparable  to  the  wavelength  of  interest  [A  ].  In  general,  because  of  the  reflections  and  refractions 
in  realistic  media,  the  wave  paths,  or  Eikonal  phase  will  be  multi-valued  in  some  regions  of  space, 
which  represent  multiple  passes  of  the  wave  front.  Then,  at  each  grid  node,  an  array  of  phase  is 
stored,  one  for  each  wave  front  passage.  This  is  easily  accomplished  using  a  counter,  which  serves 
as  an  array  index  (k].  For  each  k,  or  "trajectory",  the  recently  developed  method  -  "Wave 
Confinement"  -  (WC],  is  used  to  solve  a  modified  wave  equation.  This  method  generated  values  of 
y/k  at  grids  nodes  when  equation  (2]  is  discretized. 


drV  =  c2d2xy/  +  E  (1] 

where  E  s  a  combination  of  positive  and  negative  (stable]  dissipation.  The  purpose  of  this 
modification  is  to  generate  short,  coherent  Solitary  Waves  which  represent  wave  fronts.  When 


3  Steinhoff,  J.,  "A  New  Eulerian  Method  for  the  Computation  of  Propagating  Short  Acoustics  and 
Electromagnetic  Pulses,”  Journal  of  Computational  Physics,  Vol.  157,  2000 


equation  (1]  is  discretized,  the  solitary  waves  persist,  remaining  concentrated  over  only  2-3  grid 
cells,  unlike  conventional  numerical  schemes,  where  the  waves  dissipate.  In  this  way,  it  is 
computationally  feasible  to  solve  for  wave  propagation  over  a  large  region,  containing  many  waves. 

The  problem  described  above  is  solved  in  4  steps:  a]  Near  Field  generation  on  a  source  surface 
surrounding  known  emitters  or  targets  such  as  windmills  and  aircraft,  b]  Large  scale  far  field 
propagation,  which  computes  accurate  wave  front  propagation.  These  wave  fronts  are  much  larger 
compared  to  the  physical  waves  and  act  like  basis  functions,  whose  centroids  propagate  as  those  of 
small  scale  physical  waves.  The  Wave  Confinement  method  accurately  propagates  these  wave 
fronts  with  no  numerical  dissipation  as  nonlinear  solitary  waves,  which  retain  their  shape  and  size 
indefinitely,  c]  Computation  of  small  scale  return  signals  by  applying  the  arrival  time,  attenuation 
factor  and  propagation  vector  computed  using  b]  to  the  detailed  signal  computed  using  a].  The 
small  scale  details  are  propagated  to  the  far  field  in  this  step  by  tracing  back  the  point  of  origin 
using  the  near  field  information  saved  on  a  source  surface  around  the  scattering  objects. 

Near  Field  Scattering 

A  new  parallel  fast  Fourier  transform  and  fast  multipole  method  (FFT-FMM]  -accelerated  surface 
integral  equation  solver  has  been  developed  which  is  capable  of  analyzing  scattering  from  wind 
turbines  in  the  high-frequency  regime.  The  novelty  of  the  solver  lies  in  its  parallelization  scheme  as 
well  as  the  use  of  a  singular  value  decomposition  scheme  for  compressing  near-field  interactions. 
Flere,  the  solver  is  applied  to  the  analysis  of  scattering  from  a  wind  turbine  at  300  MHz;  the  analysis 
required  10,114,062  unknowns.  The  time  history  of  the  scattered  signal  is  computed  on  a  source 
surface  surrounding  the  near  field  of  a  scattering  object,  which  will  later  be  used  to  reconstruct  the 
signal  at  any  far  field  point.  Note:  In  the  current  project,  the  time  history  of  the  electric  field 
intensity  envelope  is  used  for  simplicity  because  the  main  objective  of  the  report  is  to  show  the 
demonstration  of  a  coupled  inner  field  near  the  wind  mill  and  the  WC  based  far  field  method. 

Far  Field  Propagation 

Since  the  medium  is  non-dispersive  and  non-dissipative,  if  we  trace  a  particular  ray  (in  our  case,  a 
short  segment  of  a  wave  front],  from  the  source  to  the  receiver  (including  refractions  and 
reflections],  the  wave  form  will  be  the  same  as  a  function  of  time,  t,  on  the  emitting  surface,  as  well 
as  at  the  receiver:  Only  the  time  will  be  translated  (by  the  travel  time]  and  the  amplitude  will  be 
multiplied  by  a  factor  that  depends  on  the  cross  sectional  area  of  the  ray  at  the  two  locations. 

Two  quantities  must  now  be  determined:  Translation  time,  and  amplitude  ratio.  A  simple  way  to 
determine  the  direction  vector  of  the  ray  that  will  reach  the  receiver  (and  hence  the  initial 
amplitude  at  the  source],  the  travel  time  and  amplitude  ratio  of  the  ray,  is  to  start  from  the  receiver 
and  propagate  backward  in  time,  taking  account  of  reflections  and  refraction.  We  employ  the  time 
reversal  invariance  property  of  the  wave  equation  to  equate  this  result  to  that  from  a  direct, 
forward  time  computation.  The  result  will  be  a  wave  at  the  source  which  is  almost  a  plane  wave, 
with  a  smoothly  varying  direction  vector,  arrival  time,  and  amplitude.  These  quantities  will  be 
computed  on  the  coarse  "outer"  grid,  which  spans  the  distance  between  the  source  and  the  receiver. 


In  general,  because  of  the  reflections  and  refractions  in  realistic  media,  the  wave  paths,  or  Eikonal 
phase  will  be  multi-valued  in  some  regions  of  space,  which  represent  multiple  passes  of  the  wave 
front.  For  this  reason  a  straightforward  finite  difference  iteration  solution  scheme  of  the  governing 
Hamiltonian  -  Jacobi  equation  is  not  feasible.  To  accommodate  these  features,  we  take  propagation 
time  as  an  extra  variable  and,  using  local  characteristics,  solve  for  the  wave  fronts  as  they  propagate 
through  these  regions,  employing  the  FDTD  scalar  wave  equation,  including  these  multiple  paths,  if 
necessary.  The  wave  fronts  that  we  describe  essentially  represent  "bundles"  of  particle  trajectories. 
The  computation  of  multiple  phase  regions  was  previously  reported  by  the  PI  in  a  closely  related 
study  of  short  acoustic  wave  propagation1. 

We  first  define  an  artificial  isotropic  wave,  y/ ,  which  is  emitted  at  the  receiver  location.  This  GF 
wave  is  numerically  confined  with  our  WC  method  and  projected  backward  in  time.  Because  of  WC, 
the  wave  remains  thin  (2-3  grid  cells]  throughout  the  propagation  to  the  distant  source  region, 
ensuring  a  close  representation  of  the  wave  front.  There,  for  long  distance  propagation,  the  GF  wave 
is  approximately  a  plane  wave.  We  assume  that  the  total  number  of  grid  cells  in  each  direction  is 
limited  to  few  hundred  for  a  feasible  computation.  For  propagation  distances  of  ~  10  km,  the  width 
of  the  captured  wave,  which  is  —  2  grid  cells  (with  WC],  can  still  be  2-3  orders  of  magnitude  larger 
than  the  desired  resolution  limit  of  the  thinnest  physical  signals.  This,  of  course,  precludes  a  direct 
numerical  solution  of  the  wave  equation,  but  does  not  preclude  a  WC  based  numerical  Greens 
function  that  retains  the  ability  to  treat  refraction  and  reflections.  WC  can  accurately  propagate  thin 
waves  since  they  follow  the  same  trajectories  as  the  much  longer  resolved  waves. 

Large  Scale 

To  summarize,  the  Greens  function  computed  using  WC  involves  computation  of  a  "GF"  wave 
"emitted"  from  the  receiver  as  i//(x,Trec  'j  and  "reverse  propagated"  to  the  physical  source.  For  an 

receiver  with  an  isotropic  receptivity,  a  unit  initial  amplitude  is  taken.  As  stated,  near  the  source, 
the  wave  amplitude,  y/ ,  can  be  multi-valued  due  to  multiple  paths  from  refraction  and  reflection. 
We  index  the  values  according  to  the  time  of  passage  in  the  (physical]  source  region.  Since  y/  is 
confined  to  only  ~2-3  grid  cells,  direct  partial  derivatives  approximations  cannot  be  used  to 
compute  derivatives  and  propagation  directions,  as  in  conventional  numerical  schemes.  Instead, 
expectation  values,  which  involve  time  integrals  as  the  wave  passes  over  each  node  (labeled  (/,/)] 


kU=Vri,j 


where  t  is  the  travel  time  of  the  reversed  wave  from  the  receiver  to  the  (/,/)  node.  An  amplitude 
ratio  ( A"e ),  which  is  the  ratio  of  the  observed  amplitude  [(y/).  ]  to  the  emitted  amplitude  is  also 

computed.  These  are  accurately  defined  due  to  the  conservation  properties  of  WC,  and  are  expected 
to  be  smoothly  varying  in  the  source  region.  The  "receiver"  times  t.J  will  also  be  smooth. 

Small  Scale 

Now,  we  project  out  the  waves  that  will  intersect  the  receiver,  k?°“rce  =  —  k,™  and  at  desired 
receiver  times,  t.e.  =  r7r“  +  T'l°“rce .  At  any  far  field  point,  the  physical  signal  is  computed  as 

E(x,trec)  =  A(x  ,trec)*E(x  ,trec) 

y  ’  )  y source ’  J  \  source’  J 

Validation 


For  validation,  an  initial  distribution  defined  by  E  =  — expy-fir2  j  on  the  source  surface  around  an 

object  is  propagated  to  the  far  field  receiver.  The  integral  variables  required  to  reconstruct  the 
physical  signal  are  computed  by  propagating  an  isotropic  wave  from  the  receiver.  This  computation 
involves  a  linearly  varying  index  of  refraction  and  could  not  be  done  using  a  Kirchhoff  method, 
which  requires  a  closed  form  Greens  function.  In  Figure  2,  the  analytical  signal  is  compared  to  the 
signal  computed  using  WC  and  they  match  quite  well. 
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Signal  at  the  observer 


Figure  2:  The  waveform  at  the  receiver  location  for  an  initial  distribution  of 


Demonstration 

Return  Signal  computation  from  a  single  wind  turbine 

The  new  method  described  in  the  above  section  is  used  to  compute  the  return  single  from  a  single 
wind  turbine  as  shown  in  Figure  3. 


Inner  Field  Computed  using  MoM 


Wind  Turbine 


Figure  3:  Problem  Setup 

Time  history  of  the  scattered  signal  computed  in  the  near  field  using  MoM  described  in  the  proposal 
is  saved  on  a  source  surface  around  wind  turbine.  The  radius  of  the  sphere  is  3R.  In  the  present 
case,  the  time  history  of  Electric  field  intensity  is  used  for  simplicity  because  the  main  objective  of 
the  report  is  show  the  demonstration  of  a  coupled  inner  field  and  far  field  method. 


Below  are  the  details  of  the  inner  field  computation  (supplied  by  Eric  Michielssen]: 

The  antennas  are  located  few  kilometers  away  from  the  wind  farms  and  so  the  EM  wave  that 
reaches  any  wind  turbine  is  planar.  So,  the  wind  turbine  is  illuminated  by  a  plane-wave 
(propagating  along  y  direction].  E-field  with  unit  magnitude  is  polarized  along  z  direction.  The 
blades  of  wind  turbine  lie  on  yz  plane  and  the  ground  is  assumed  to  be  located  on  xy  plane.  The 
near  field  points  are  selected  on  a  sphere  with  the  radius  of  125m  (3  times  blade's  length].  Center  of 
sphere  is  positioned  at  (0,  0,  0],  50  samples  are  selected  along  theta  direction  and  100  samples  are 
selected  along  phi  direction.  Intervals  for  theta  and  phi  angles  are  assigned  as  [0  180],  [-90  90],  (in 
degrees],  respectively.  So  totally,  there  exist  5000  near  field  points  on  half-sphere.  The  period  of 
blade  rotation  is  set  to  3  seconds.  The  envelopes  of  sums  of  squares  of  E-field  components  at  near 
field  points  are  computed  at  1441  points  from  t=0  to  t=3  seconds  and  stored  on  the  source  surface. 
In  Figure  4,  a  snapshot  of  the  Electric  field  intensity  is  shown  on  a  sphere  of  radius  3  rotor  radii 
with  wind  turbine  at  the  center. 


Plane  wave 


Figure  4:  Near  field  computation 

Nonlinear  solitary  wave  propagation,  taking  into  account  reflection  and  refraction,  can  provide 
information  about  the  region  on  the  inner  field  that  accounts  to  the  signal  at  the  receiver,  which  will 
then  be  used  to  compute  the  signal.  For  multiple  arrivals,  point  of  origin  and  return  signals 
corresponding  to  each  arrival  are  computed.  For  the  above  case,  the  return  signal  at  a  receiver 


placed  ~700m  away.  The  computed  first  arrival  (incident]  and  second  arrival  (reflected]  signals  are 
plotted  against  the  exact  signals  in  Figure  5  and  Figure  6  respectively  and  they  seem  to  match  well. 
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Figure  5:  First  Arrival  for  One  Windmiliat  ~700m  distance 
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Figure  6:  Second  Arrival  for  One  Windmill 

Return  Signal  From  Multiple  Wind  Turbines  With  and  Without  Aircraft 


Another  demonstration  to  compute  the  return  signal  from  a  notional  wind  farm  with  and  without 
aircraft  is  described  below.  The  locations  of  receiver,  wind  farm  and  aircraft  are  shown  in  Figure  7 
and  a  schematic  of  the  problem  description  is  shown  in  Figure  8. 
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Figure  7:  Locations  of  emitters 
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Figure  8:  Location  of  wind  turbines  and  aircraft  on  xy  plane  at  height,  z  -  0 


The  computation  is  performed  with  and  without  aircraft.  The  combined  signal  due  to  multiple 
arrivals  is  shown  in  Figure  9  and  it  can  be  seen  that  the  computed  signal  compares  well  with  the 
exact  signal. 
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Figure  9: 10  Windmill  Farm  Radar  Envelope  over  1  Windmill  Period 
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Figure  10:10  Windmill  Farm  Radar  Envelope  over  1  Windmill  Period,  with  Airplane 
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Figure  11:  Wind  farm  Envelope  with  and  without  Airplane 


